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Abstract 


Two approaches to determining the ocean sound speed profile using measured acoustic modal 
eigenvalues are examined. Both methods use measured eigenvalues and mode dependent 
assumed values of the WKB phase integral as input data and use the WKB phase integral 
as a starting point for relating the index of refraction to depth. Inversion method one is 
restricted to monotonic or symmetric sound speed profiles and requires a measurement of the 
sound speed at one depth to convert the index of refraction profile to a sound speed profile. 
Inversion method two assumes that the sound speed at the surface and the minimum sound 
speed in the profile are known and is applicable to monotonic profiles and to general single 
duct sound speed profiles. For asymmetric profiles, inversion method two gives the depth 
difference between two points of equal sound speed in the portion of the profile having two 
turning points, and in the remainder of the profile it gives sound speed versus depth directly. 
A numerical implementation of the methods is demonstrated using idealized ocean sound 
speed profiles numerical experiments used to test the performance of the inversions using 
noisy data. The two methods are used to determine the sediment sound speed profiles in 
two shallow water waveguide models, and inversion method one is used to find the sediment 
sound speed profile using data from an experiment performed in the Gulf of Mexico. 
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Chapter 1 


Introduction 


As sound propagates through a medium, the sound field is affected by the medium through 
which the sound is passing, and by the nature and physical properties of the medium’s 
boundaries. This relationship between the sound field and the medium may be used to predict 
the properties of the sound field as the sound propagates through a medium with known 
properties (the forward problem) or to infer the properties of a medium from measurements 
of the sound field after propagation through the medium (the inverse problem). Using sound 
we may thus able to measure the properties of media which might otherwise be difficult 
to probe. In this thesis we present a nonperturbative inverse method based on the WKB 
approximation to the sound field which uses the values of the modal eigenvalues as input 


data. 


1.1 Background 


In ocean acoustic tomography, linear inverse techniques are applied to measurements of per- 
turbations in the sound travel time which result from sound speed variations along the signal 
path. These sound speed variations are then related to the oceanic processes causing the 
sound speed variations. Acoustic tomography techniques have been used to study ocean fea- 
tures such as mesoscale eddies [12], internal waves [17], and barotropic motions [9]. While 
these processes may be studied in other ways, tomography has the advantage of being able 


to provide frequent measurements of the ocean properties over large areas. Equivalent mea- 
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surements made by other means would be prohibitively expensive. In these tomographic 
studies, measurements are made using sound travelling over paths which do not involve in- 
teraction with either the ocean surface or bottom, however; sound which does interact with 
the boundaries can be used to study the boundaries. James Miller in his ScD thesis [34], 
for example, uses tomographic techniques to obtain estimates of the sea surface spectra, and 
acoustic signals generated in the water column have been used extensively in marine seismol- 
ogy to study crustal structure and properties [16]. For studying sound propagation in the 
ocean, the properties of the ocean bottom are needed on a much finer scale than provided 
by marine seismological measurements, and the region of interest, rather than extending to 
several kilometers depth, is limited to the top few hundred meters of the sediments at most 
frequencies of interest. Frisk et al. [19] use amplitude versus range data obtained using a 
deep-towed pulsed CW source and two receivers moored near the bottom to infer geoacoustic 
models of the upper sediment layers in deep water. The model is derived from the measured 
data using an iterative forward modelling technique. This method relies on time gating to 
separate the bottom reflected signal from the surface reflected signal. In shallow water, time 
gating method is not practical because the time differences between signals reflected from 
the surface and signals reflected from the bottom are too small. Instead, a horizontal array is 
used to measure the steady state sound field which is then decomposed into its eigenfunctions. 
Using the eigenvalues of the propagating sound field, Rajan et al. [37] obtain the sediment 
properties using a perturbative approach. 

Rather than using a perturbative technique, we use a WKB phase integral based inverse 
method to obtain the sound speed profile using as input data the eigenvalues of the sound field 
in the ocean acoustic waveguide. The method may be used to find the sound speed profile in 
the water column and in the upper layers of sediments. Although the sediments have enough 
rigidity to transmit shear waves, and we expect compressional to shear wave conversion to take 
place at layer interfaces, Fryer [23] indicates that effects of compressional to shear conversion 
resulting from shear speed gradients within the sediment column is small above 20 Hz. Since 
we are interested in the upper sediment layers where shear wave speed and shear wave speed 
gradients are small, we will treat the sediments as a fluid extension of the water column and 


neglect all shear effects. This will enable us to determine the compressional sound speed 
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profile in the sediments in that portion of the sediments affecting sound propagation at the 
frequencies of interest (the standard frequency used in the examples will be 220 Hz). These 
methods will not be able to distinguish the low velocity zone that commonly occurs at the 
water /sediment interface for bottom materials such as pelagic clay [25]. This inability to find 
low velocity regions results because the inversions are determining ray turning depths and 
rely on sound energy being turned by an increasing sound speed within a region in order to 
find the sound speed/depth relationship in that region. The sound does not turn within a 
low velocity zone. Using an approach suggested by Vassell to study the index of refraction 
in optical waveguides [43], we relate the normal mode eigenvalues to the sound speed profile 


via the WKB phase integral 
22 1/2 
ko | [n?(z) — sin”, | ; ee ie) 
21 


Here ko is the wavenumber given by w/co where w is the frequency and co is the reference 
sound speed, @, is the local propagation angle corresponding to the nth mode, and z, and 
zq are the turning depths for the mode. The phase integral will be treated as a function of 
the new variable € = sin*6,. The value of this integral may be estimated based on the mode 
number and the type of boundary interactions experienced by the mode. For example, the 
phase integral for the third mode will assume a different value in an environment giving two 
turning points than it will in an environment where the mode interacts with the surface and 
has only a lower turning point. 

There are two methods for obtaining the index of refraction as a function of depth. The 
first method (inversion method 1) follows Vassell’s treatment and can be used for sound speed 
profiles that increase monotonically with depth or are symmetric about a depth of minimum 
sound speed. In this method, the WKB phase integral equation is differentiated to arrive at 
an Abel integral equation which can be used to relate depth to index of refraction. From 
knowledge of the sound speed at one depth, we can determine the sound speed profile. If the 
sound speed profile is symmetric rather than monotonic, the inversion gives the sound speed 
profile below the depth of symmetry, and the symmetry of the profile is used to generate the 
upper portion of the profile. 

The second method (inversion method 2) is based on a method used in quantum me- 


chanical scattering for determining molecular potentials and is referred to in the quantum 
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mechanical context as the JWKB or semi-classical approximation [10],[44]. Now the WKB 


phase integral is transformed to the form 


He= ["[w(@)-g de (1.2) 


and J(€) (referred to as the inclusion) is differentiated to produce the excursion which is the 
difference in depth between two points of the sound speed profile having the same sound 
speed. In general there is insufficient information to translate this to a sound speed profile; 
however, if the profile is monotonic or symmetric, this method gives the sound speed profile 
directly. For an asymmetric profile in which some modes have two turning points and some 
have one turning point, inversion method 2 will give a sound speed profile for the range of 
sound speeds having a single turning depth since the surface is essentially the upper turning 
depth. 

Munk and Wunsch [35] derive an Abel transform based inversion method along the lines 
of our inversion method 1; however, they deal strictly with applications to the water column, 
and they approach the asymmetric case in a different fashion. Much of the previous work in 
inverse methods for finding the sound speed profile in the ocean and the sediments has used 
perturbation methods largely because the techniques of linear perturbation theory provide a 
means to compute estimates of errors and depth resolution [33]. With the methods discussed 
here, we can only give qualitative arguments concerning resolving power and numerical ex- 
amples illustrating accuracy of the methods and their performance in a noisy environment. 
The main advantage with our methods is that they do not require an initial estimate of the 


sound speed profile. 


1.2 Overview 


Chapter 2 contains a review of the essentials of normal mode theory, ray theory, and WKB 
theory in the context of underwater acoustics. Since the inversion methods presented here 
are based on the normal mode eigenvalues as input data and use the WKB phase integral 
as a Starting point for relating the index of refraction to the eigenvalues, an understanding 
of normal modes and the WKB approximation is important. Ray theory provides a means 


for building a conceptual picture of the inversion process through the ray theory view of the 
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turning point. 

Chapter 3 contains the bulk of the thesis. We start by grouping normal modes according 
to the boundary interactions and the number of turning points (we do not allow for multiple 
duct profiles), and then evaluating the WKB phase integral for each of the four types of 
mode. These phase integral values and the associated eigenvalues provide the input data used 
in generating a profile dependent functional relationship for the phase integral. The Abel 
integral equation based inversion relations for the two inversion methods are then derived. 
Following a brief description of the spline methods used in the computations, we discuss 
the numerical implementation of the two inversion methods. Both inversion methods are 
applied to the n?(z) linear profile, the parabolic profile, and a bilinear profile to illustrate 
their capabilities under ideal conditions. Using the n*(z) profile, we then demonstrate the 
performance of the inversion methods using data contaminated with random noise. 

In chapter 4 we apply inversion method 1 to the problem of determining the sediment 
sound speed profile in a shallow water waveguide by treating the sediment as a fluid extension 
of the water column. Two synthetic models are used to test the inversions: (1) a waveguide 
with a terrigenous bottom, and (2) a waveguide with a fine sand bottom. We then apply the 
inversion method to data obtained in an experiment performed in the Gulf of Mexico. 


Chapter 5 summarizes the results of the work and presents ideas for future work. 
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Chapter 2 


Normal Mode Propagation and the 
WKB Approximation 


2.1 Modal Propagation 


The input data for the inverse method considered in this work are the measured horizontal 
wavenumbers or eigenvalues for the modal description of sound propagation in the medium, 
and a brief review of normal mode propagation is given here to emphasize the normal mode 
ideas important to the inversion. Normal mode theory provides a ‘full wave’ solution to 
the wave equation in the ocean waveguide which takes into account the properties of the 
medium in the waveguide, frequency effects, and the effects of the properties of the waveguide 
boundaries. The degree to which the solution to the wave equation provided by modal 
theory realistically describes acoustic propagation in the ocean depends on the accuracy to 
which the medium properties and boundary conditions are modeled. We will use models of 
varying sophistication to illustrate various aspects of the normal mode description of sound 
propagation pertinent to the inverse method. 

We model the ocean as a waveguide consisting of a water column bounded above by 
the air-ocean interface where the acoustic pressure is zero, and below by a sediment bottom 
which will be considered a semi-infinite half-space. The sound speed is c(z) and the density 
p are taken to be piecewise constant throughout most of the work. The following simplifying 


assumptions are used throughout: 
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1. The source is a monochromatic point source. 
2. Cylindrical symmetry about the vertical axis through the source. 


3. The medium properties vary only in depth, i.e. the problem is range independent. 


Included here is the assumption that the water depth is constant. 
4. No scattering from rough surfaces or inhomogeneities in the medium. 


5. The bottom is a fluid which may effectively be treated as an extension of the water 


column. 


With the assumption of cylindrical symmetry, the spatial part of the acoustic pressure 
p(T, z, 20), due to a point source located at r = 0 and z = zo with harmonic time dependence 
exp(—iwt), satisfies the inhomogeneous Helmholtz equation: 

10 0 0? d(r 
Eo (=) + 72 4. (2) p(T, 2,20) = —2 2) 6(z — 20) (au) 
where w is the frequency and k(z) = w/c(z). 

To complete the problem, appropriate boundary conditions must also be specified. At 
any interface where there is a sudden change in the medium properties (the water-sediment 
interface for example) the boundary conditions of continuity of pressure and continuity of 
normal particle velocity must be satisfied. In the range equation the radiation condition that 


there are no waves propagating inward from infinity will be used. 


2.1.1 The Rigid Bottom 


We start by considering the simplest model of the ocean waveguide consisting of a layer of 
water with sound speed c(z), which may be a function of depth, and constant density, po. 
The layer is bounded above by a pressure release-surface (the air-water interface) and below 
by a rigid bottom. The symmetry and range independence of the problem lead naturally to 


use of separation of variables in solving the problem. Letting 


P(7, 2,20) = $(z)¥(r), i22)) 
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and substituting this into the homogeneous Helmholtz equation leads to the separated equa- 








tions 
d* $(z) 2 
ge TN _ Zs 
qe (1 OZ) = 0 | (2.3) 
and 
d*y(r) | 1dy(r) 2 
dp sts aap +k v(r) = 0 (2.4) 
where —k? is a separation constant, and 
~? = k?(z) - ae (2.5) 


The total acoustic pressure vanishes everywhere on the pressure-release surface and a 
sound wave incident on the surface is completely reflected with a plane wave undergoing a 7 
phase shift. The rigid boundary condition is characterized by the condition that the normal 
particle velocity vanishes on the surface resulting in complete reflection. For an incident plane 
wave there will be no phase shift upon reflection from the rigid boundary. In this model, the 


depth function $(z) satisfies the boundary conditions 


(0) = 0 (2.6) 
for the pressure-release surface and 
dp(h) _ 4 (2.7) 
dz 


for the rigid bottom. 

Equation (2.3) only has solutions for certain values of the parameter y (the eigenvalues). 
Because this is a Sturm-Liouville problem, the set of eigenfunctions ¢,;(z), 7 = 1,2,3,..., 
forms a complete orthonormal set of square integrable functions on the intervalO < z<h 
implying that [45] 

footed) 0 4 Bs; 
0 PO a) 
where the 1/9 weighting is included in analogy to the results in problems where the density 
is a function of depth. 
Up to this point we have allowed the sound speed to be any function of the depth z. For 


definiteness in illustrating the solution of the problem, we now specify c(z) = constant. The 
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boundary conditions lead to the equation 


P(g Mes 
a2? 


for the vertical wavenumber y,. Using a density pp = 1 and normalizing gives for the 
eigenfunctions 


Onze sin |(k? ne ee] (2.9) 


where the horizontal wavenumbers k,, are given by 


k? = k? — (et (2.10) 


Expanding the solution of equation (2.1) in terms of the eigenfunctions 
Le @) 
P(r, 2,20) = >, bn(z)n(r), (Cail) 
n=] 
substituting expansion into equation (2.1), and using the orthonormality of the depth func- 


tions gives 
d?pn(7) Le ldyp,(r) 
dr? r dr 


which is a Bessel equation of order zero for which the solution is a linear combination of Hankel 


+ k? ¥alr) = -=6(r)bn(20) (Zale), 


functions. We use the radiation condition that there are no waves propagating inward from 


infinity to eliminate the term involving Ro), and the solution to the radial equation is 
Wall) = ibn (zo) HO )( ke omy (2503) 
The solution to equation (2.1) is then expressed as a sum of the normal modes 
Dye) (1) 
P(r, Z,20) = 7 yy SIM(¥_ 2) SIN J,,20)) ge eee): (2.14) 


At distances of more than a few wavelengths from the source, we can use the asymptotic 


approximation to the Hankel function [2] 


HH" (k,r) ~ Ja expli(k,r — 7/4)] (2.15) 


The normal mode solution of equation ( 2.1) given by equation ( 2.14) describes the 





sound field in terms of a sum of travelling cylindrical waves propagating outward from the 
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origin. Each of these modes is a standing wave in the vertical direction with a distinct depth 
distribution of pressure described by the eigenfunction ¢(z) and the vertical wavenumber 7,. 
Although the summation of equation (2.14) is infinite, for a sufficiently large mode number N, 
the horizontal wavenumber will become imaginary, and the solution will rapidly decay with 
range. For a receiver located sufficiently far from the source, the sound field will be adequately 
described by a finite sum of propagating modes. It is these horizontal wavenumbers from 
this description of the sound field that will provide the input information for the inversion 


method described in the next chapter. 


2.1.2 Normal modes in a waveguide with a penetrable bottom 


If the bottom of the waveguide is not a rigid reflector, some of the sound energy will penetrate 
into the bottom, and the solution to the Helmholtz equation must be found both in the water 
and in the bottom. The approach taken here is different from the previous section in that 
the Hankel transform is used to obtain the solution; however, the essential ideas related to 
the propagating modes remain the same. 

Starting once again with the inhomogeneous Helmholtz equation (equation ( 2.1)), we 


introduce the zero order Hankel transform pair 
CO 
P(r, 2,20) = / D(kr, 2, 20) Jo(krr )krdr (2.16) 
0 


(hr, 2,20) = f p(r, z, 20) Jo(krr)rdr a) 
0 


where Jo(k,r) is the Bessel function of order zero. 

Multiplying equation (2.1) by Jo(k,r)rdr and integrating from 0 to oo (see [3] for details 
of the integration), we find that the transform variable (or depth dependent Green’s function) 
must satisfy [20] 

He 


(= 4. k7(z) — i) p(k, Z, 2%) = —26(z — 20). (2:18) 


Let p, and pz be linearly independent solutions of equation (2.18) chosen such that p, satisfies 
the surface boundary condition, and p> satisfies the bottom boundary condition. Then the 


transform function p(k,, z) is given by [20] 


oy ky, D k, 
Blk, 220) = PC AEE) oc scx (2.19) 


19 


= =—27 Kes D Ky, 2 
P(k,,z, 20) = EO El, ae (2.20) 


where W(zg) is the Wronskian given by 
W (20) = po(kr, 20) By (kr, zo) _ po(k,, z0)p1(kr, Zo) (2.21) 


where the prime indicates differentiation with respect to z. 


Next the Bessel function is expressed in terms of Hankel functions as 
1 
Jo(krr) = 5 HS” (ker) + HS” (kyr)| (2.22) 


and use the relation 


H)(e-'* kyr) = —~HY (k, 7) (2.23) 


to allow the range of integration for the transform integrals to be extended over the whole 
real k, axis from —oo to oo when evaluating the transform integral for p(r, z). The solution is 
then found by contour integration methods [3], and may be expressed as the sum of a discrete 


and a continuous portion given by [20] 
p(r,z) = it) bn(z0)On(z) Hb” (kent) + I(r) (2.24) 


where the eigenfunctions and eigenvalues satisfy 


(+ +k*(z)- o gn(z) = 0 (2.25) 


along with the surface and bottom boundary conditions. 

The part of the solution representing the continuum (J(r) or the branchline contribution) 
is rapidly attenuated with range and will be neglected, since we are only concerned with 
measurements made more than a few wavelengths from the source. We note that the dis- 
crete portion of the solution is of the same form as for the rigid bottom case, and that the 
eigenvalues have the same interpretation, although they will take on different values since 
the sound is propagating in both media, and the characteristic equation will be different.In 
general the lowest medium will be taken as a half space extending to infinity and a radiation 


condition will be imposed resulting in a solution that decays with depth in the halfspace. 
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2.1.3 Analytically solvable profiles 


We next consider two specific sound speed profiles in order to illustrate the concept of turning 
depth. A turning depth is a depth at which the vertical wavenumber vanishes, leading to 
a fundamental change in the behaviour of the solution to the Helmholtz equation. In the 
region where the vertical wavenumber is greater than zero, the solution is oscillatory, while 
in the region for which the vertical wavenumber is less than zero, the solution will decay 
with depth. Depending on the sound speed structure, these turning points may be located 
in either the water column or in the sediment. The sound speed profiles we will consider are: 
(1) the square of the index of refraction a linear function of depth, and (2) the square of the 
index of refraction parabolic. In both these cases solution of the depth equation gives closed 
form expressions for the eigenvalues. Additionally, we will later find that these profiles may 
be inverted analytically by the methods presented here, and consequently they will be the 


standard examples for illustrating results and verifying the numerical methods. 


The n*(z) linear waveguide 


For simplicity we consider a semi-infinite ocean with a pressure-release surface and an index 


of refraction which satisfies the relation 
n*(z) = 1—az. (2.26) 


In lieu of a bottom boundary condition, we apply a radiation in depth requiring that the 
solution remain finite at infinity. 
Proceeding as for the penetrable bottom case, the homogeneous part of equation (2.18) 


becomes 
d* p(k, 97% zo) 


wate + [h5(1 - 2) — ke] Br, 2,20) = 0 (2.27) 


Making the change of variables [8] 
H =(ak2)-*3. to = H?(k?2 — 2); t=to+2/H, (2.28) 


equation ( 2.27) is transformed into the Airy equation 





THY) _ 1501) (2.29) 


a 


which has as solutions the Airy functions Ai(t) and Bz(t) [2]. (Here functional dependence 
other than dependence on t has been surpressed.) We eliminate Bi(t) in order to satisfy the 
radiation condition, and the eigenvalues are found using the pressure release condition which 
requires Ai(to) = 0 from which we see that to = —Yn where the yn are the zeros of the Airy 


function. With t, = z/H — yn and ton = 2o/H — y, the sound field is 


it —~ Ai(t, Ailton Ho (ker 
P(r2) = Fr | ao = -_ 


The eigenvalues are given by 
ke hee ie 231) 
The parabolic profile 


In considering the parabolic profile, we take the ocean as infinite in the vertical direction 
so that there are no boundaries, and the radiation condition is used to give a solution that 


remains finite at too. The index of refraction satisfies 
n*(z) = 1—-—a*2’, (2232) 


and the homogeneous part of equation (2.18) is 


daphne zy z0) 


de a [aC are k?| (co) 0 (2.33) 


The solution to this are given in terms of Hermite polynomials as [42] 
D Reg) (2"n!)~"/2 (ako /m)*/4e7 2%?" H (Jako Zi) (2.34) 
where the H,,(-) are the Hermite polynomials. The eigenvalue condition gives 


kn = ko (1 —a(2n+1)/ko)/?7; =n =0,1,2,.... (2.35) 


2.2 Rays 


Ray theory is a geometric solution of the wave equation which is correct for high frequencies 
and provides a simple physical description of the solution in terms of the paths along which 


the acoustic energy is refracted. The ray approximation is widely used for studying sound 


ze 


propagation in underwater acoustics, and provides a convenient method for visualizing the 
propagation paths of the sound energy and the depth of penetration of the sound energy into 
the bottom. We will briefly illustrate the connection between the normal mode view and the 
ray picture to the extent that it helps in understanding the inversion technique. A number 
of references are available which provide detailed treatments of ray theory in underwater 


acoustics [7],[6]. 


2.2.1 The Equations of Ray Theory 


The equations for ray acoustics [7] are obtained from the Helmholtz equation in vector form 


using pressure as the variable 


V7p+ k*(z)p = 0, (2.36) 


where Z is the position vector and k(z) = w/c(z), by taking the solution to be 
p(t) = A(z) exp(tkyS(z)). (2:37) 


Here ko = w/ Co where Co is a reference sound speed such as the minimum sound speed in the 
waveguide (or the sound speed at the source), A(Z) is the amplitude of the sound waves, and 
koS(z) is the phase, with the function S(z) referred to as the eikonal. 


Substituting the assumed solution ( equation ( 2.37)) into the Helmholtz equation gives 
V?A + iko(2V.A- VS + AV?S) + K3A [n?(z) — (VS)"] = 0. (2.38) 


The equations of ray theory are obtained from equation ( 2.38) in the limit as kg — oo by 
equating the real and imaginary terms to 0, and neglecting the V?A compared to the real 
term containing kA (that is we assume V*A/A < k@). The resulting equations are the 


eikonal equation 


7S) ear) (2.39) 


and the transport equation 


2VA-VS+AVS =0. (2.40) 


The eikonal equation defines the geometry of the rays which are lines perpendicular to the 
surfaces of constant phase, that is S(z) = constant . These ray trajectories can be computed 


to determine the path of the sound energy. 
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2.2.2 The Connection Between Rays and Modes 


C.T. Tindle and K.M. Guthrie [24] [41] have illustrated the connection between normal mode 
theory and the ray theory approximation. Considering a deep water sound speed profile with 
a minimum sound speed at some depth below the surface, then for a given mode number n, 
a ray with the same turning point depths as mode n can be defined by 


Co Ww 
— = — Ay 
sin@o =i, ( ) 


where the angle 99 is the angle of incidence for the ray as it crosses the channel axis (i.e. the 
sound speed minimum), and k,, is the horizontal wavenumber for mode n. The ray effects 
are duplicated by summing over a group of modes, and each group will manifest itself as an 
energy pulse travelling along the ray path which corresponds to the ray defined by the angle 
$9 for the central mode in the group of modes. Thus we can construct a ray corresponding 
to the group of modes; however, this ray will not necessarily be the same as the physical 
ray path along which the sound energy actually travels between the source and receiver. By 


using this view of rays as the result of modal interference, we can construct a picture of the 


path followed by the sound, and its turning depths within the framework of normal modes. 


2.3 The WKB Approximation 


WKB theory, like ray theory, provides an approximate solution to the wave equation in 
the high frequency regime; however, the WKB method takes frequency information into 
account in the amplitude function as well as in the phase, and it accounts for frequency 
dependent effects such as diffraction which are ignored in standard ray theory. In general, the 
WKB method is a means for finding an approximate solution to a linear ordinary differential 
equation in which the highest order derivative is multiplied by a small parameter [5] and its 
use is not limited to underwater acoustics. In the problem at hand, the differential equation 
of interest is the equation for the depth dependent Green’s function (equation ( 2.18)) 


d*D( kr, zy zo) 


d*z a ko[n?(z) a k? /k5]b( kr z, 20) = 
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where ko = w/co, co is the sound speed minimum, and n(z) is the index of refraction. For 


later convenience we make the following change in variables 
K(z)=n%(z),  €= RRR (2.42) 
and 
Q?(z) = w(z) -& (2.43) 
Equation (2.3) then becomes 


FAK 2120) 5 6807(<)6(ke, 2,20) = 0. (2.44) 


This equation is solved by expanding 7 as (22] 


D(kr, Zz, zo) = erp kof dz ey v(erhs*| (2.45) 
Ay v=0 


where the lower limit z is some conveniently chosen constant depth. Substituting this ex- 
pansion into equation (2.44 ) and setting the coefficients of successive powers of ko equal to 


zero, we obtam 


a. (2.46) 

dy,— y 

— = - » hear ese | 7 ae (2.47) 
p=0 


from which the series in the exponent of equation (2.45) may be determined. This series will, 
in general, be asymptotic and not convergent. Retaining the first three terms and using a 


convenient normalization gives the second order WKB approximation 
a 1 
Ale, 2,20) = (hoQ(2))""? exp | f” (1+ =~) boQ(=)de (2.48) 
20 0 
where 
3/2 a 7 
€o = (Q(z)ko) >” qe? (koQ(z))-”?. 
The steps leading to equation ( 2.48) are justified if 


| yoko? |<] mko' |< | yo | (2.49) 


OT 


L601 <| & (bo@(2)) | «2 (2.50) 
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| / OOO ee (2.51) 


then we can make the physical optics approximation [5] of retaining only the first two terms, 


and the two linearly independent solutions of equation (2.3) are 
Blk, 2,29) = ky! («(z) — 6)! exp E iko / («(z) — €)'/? dz}. (2.52) 
20 


The WKB method has provided an approximate solution to equation (2.18) which takes 
into account sound speed profile and frequency information, and which allows for frequency 
dependent effects such as diffraction. These WKB modal solutions can then be used in the 
modal sum (equation (2.24)) to approximate the full wave solution. It is well known that the 
WKB solution fails at the turning points i.e. the points at which Q(z) = (K(z) — €) =0. As 
noted earlier, the behaviour of the solution differs on the two sides of the turning point with 
oscillatory behaviour above the turning point and an exponentially decaying behaviour below 
the turning point. There is an extensive literature concerning the problem of connecting the 
solutions through the turning point (see [5]); however, we are only interested in the phase 
memory of the solution i.e. the integral in the exponential for the region above the turning 


point, and not in generating approximations to the sound field. 
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Chapter 3 


The Inversion 


We approach the inversion by first grouping modes according to the boundary interactions 
they undergo and the number of turning points. Based on these factors the WKB phase inte- 
gral is evaluated for each type of mode resulting in an expression which depends only on the 
mode number in the case of modes which do not interact with the bottom, and on mode num- 
ber and the phase of the bottom reflection coefficient for bottom interacting modes. Treating 
these expressions as functions of the variable £, we proceed, for each inversion method, to 
obtain an Abel integral equation which enables us to extract a relationship between the index 
of refraction and either the depth (inversion method 1) or excursion (inversion method 2). 
Both inversion methods are applied to three prototypal sound speeed profiles to demonstrate 
their performance under ideal conditions. The performance of the inversions is then tested 


using one profile (the n*(z) profile ) with random errors added to the eigenvalues. 


3.1 Evaluating the Phase Integral 


There are four basic types of normal modes depending on the region where the normal modes 
are concentrated i.e. the region where \/k — € is real [7] (sound speed profiles with multiple 


sound ducts are excluded from this discussion): 


1. 0<2z< 2,,,.,- Lhe region is bounded above by the pressure release surface and below 


by a turning point. 


at 


2. Ztupper < 2 < Ztiower: Lhe region is bounded above by a turning point at z,,.., and 
below by a turning point at z,,...- On the ocean surface and bottom, the sound field 


is exponentially small. This case would include, for example, SOFAR propagation. 


3. 0< z< H. The region is bounded above by the ocean surface and below by the ocean 


bottom. 


- Ztupper < 2 < H. The region is bounded above by a turning point and below by the 


ocean bottom. 


Because we are interested in the properties of the sediments as well as the water column, we 
include in the categories having lower turning points, those instances in which the sound is 
not totally reflected at the interface, but turns within the sediment layers. 

For convenience, let 24, = 21,2. 2t2 = Ztupper? and y = (K — £)}/2, The WKB solution to 


the depth equation for the transform variable p(z) is 


By = Creep (ike e dz) + Cyexrp (~ike ydz)) 


Bee zy (Sl) 


C3exp (to | + | ts) , 2h es (3.2) 


Zt, 


eae l 
P(z) = M71 


The two terms in equation (3.1) represent upgoing and downgoing waves in the sound channel 





while equation ( 3.2) represents an exponentially attenuated wave below the turning depth. 


The constants are related by (8] 
Ci =Cz3exrp(ix/4); C2 = Cz3erp(—i7/4) 


with the 7/4 factor accounting for the phase associated with turning at the turning depth. 
Our interest is above the turning depth where the solution is oscillatory, and that solution 
can be written 
ya 2C3 ae 
p(z) = —cos (io f 4 dz — x/4) (323) 
( 7 ‘ 
from which the phase integral will be evaluated using the boundary conditions. We are not 


interested in the value of the constants C3, and will deal only with the phase integral. 
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For the first mode type (case 1), the pressure release condition at the ocean surface 


requires that p(0) = 0. Setting the cosine term equal to zero requires 
Zt, 
ko | Yn dz = (n—1/4)7; alae (bec) 
0 


For the modes with two turning points (case 2), we assume that the turning points are far 
enough away from the boundaries that interactions with the boundaries are not of concern. 


Expanding the cosine into the sum of exponentials 


Zt 1 ee 
COS (xo | ‘Yn dz — in/4) = 5e2p (ito f Yn dz — in/4) 
z cS 


= =eap (—iko 7 Yn dz + in/4) ; (3.5) 


The first term is a wave propagating in the positive z direction (down), while the second term 
is a wave propagating in the negative z direction (up) i.e. a wave which was reflected at the 
upper turning point. At the turning depth, the wave is reflected with a 7/2 phase change. 


Taking the ratio of the two terms and setting the ratio equal to exp [iz(2m + 1)] gives 
ko | a, dz = (n —1/2)7; ie Aa. eee (3.6) 
2t9 


with m = (n — 1) used to provide indices starting at 1. 

For the case 3 we consider a mode interacting with the pressure- release surface above 
and totally reflected at the bottom. Expanding the field into up and downgoing waves again 
gives 


ple) = Aeep (ito [ Yn iz) + Bezp (—ito | aye dz) : (3.7) 
0 0 


Evaluating the terms at the two boundaries in terms of the surface and bottom reflection 





coefficients 
p(0) | _ A 
nye” eet: 3.8 
0) 1 B = 
and 
p(H)t _ B 7 [ 
R= =— —2 az A 3.9 
BS erp Ala le DO (3.9) 
Eliminating the ratio A/B leaves 
H 
R,R,exp (2i / oP ts) = 1 (3.10) 
0 
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Letting R, = exp(12®,) and R, = exp(i2®,), gives 


ee 
%, +o, + / Vaz — (ll es ie RO ne (ol) 
0 
For a pressure release surface ®, = —71/2 and 
H 
/ Fade ae 1/2 = (3.12) 
0 


where ®, is the phase of the reflection coefficient for total internal reflection at the angle of 
incidence for the given mode. In the case of the Pekeris waveguide ( an isovelocity water 


column over an isovelocity bottom half-space) this is [11] 


t= nrc taro (3213) 
P1C1 COs A 
b, = (2) in? = (3.14) 
co 


where the p are the densities, c are the sound speeds, and @ is the angle of incidence in the 
upper medium. The subscript 0 refers to the upper medium, and the subscript 1 refers to the 
lower medium. 

In case 4 the wave reflected at the upper turning depth lags in phase by 7/2 behind the 
incident wave [7], and ®, = —7/4 giving for the phase integral 


H 
ko | Yn dz = (n — 3/4)x — 9); i as (Salk) 


3.2 Derivation of the Inversion Relations 


We designate the WKB phase integral divided by the wavenumber ko as F(£). The variable 
€ is physically equal to sin? 6 where @ is the angle of incidence of the ray corresponding to 
the mode having the eigenvalue kp,/f. Viewing the modes as a consequence of constructive 
interference of rays, we treat € as a continuous variable with the measured values of the 
horizontal wavenumber and the WKB phase integral defining values of the function F(£) at 
selected points. An Abel integral equation which can be used to relate depth to the index 
of refraction (or sound speed depending on the normalization used) may be obtained by two 
methods. The first approach is applicable to sound speed profiles which are monotonically 


increasing with depth or symmetric with respect to the sound speed minimum. The second 
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approach allows determination only of the distance between points on the two branches of 
the sound speed profile having the same sound speed when the profile is asymmetric. If some 
of the modes for the asymmetric profile have two turning points and some have one turning 
point, then for sound speed range for which there are two turning points the result will be as 
described, but for the sound speed range with one turning point, the result will be the sound 
speed profile. For a sound speed profile which is monotonically increasing with depth or is 
symmetric about the sound speed minimum, the second method, like the first, provides depth 
as a function of index of refraction . We will refer to the first method which is applicable 
to monotonic and symmetric profiles as inversion method 1, and we will refer to the second, 


more general, method as inversion method 2. 


3.2.1 Derivation of inversion method 1 


Considering the case of a profile which is monotonically increasing with depth, the phase 


integral of the WKB solution to the wave equation (assuming constant density) is written as 


[PO - Galo)" are Be ra z aie 


Here ko is w/co, co is the value of the sound speed at the surface, and k,, is the horizontal 
wavenumber for the nth mode. The equation has been normalized by dividing through by ko, 
and the quantity which will ultimately result from the inversion is the index of refraction n(z). 
Here z; is used to denote the single turning depth. If the sound speed profile is symmetric, 
the phase integral is taken between upper and lower turning points, and the symmetry of the 
problem allows the integral to be split into two equal pieces with the resulting factor of two 
taken to the right-hand side of the equation. 


With the change of variables from the previous chapter 
K(z)= (2);  € = (ra /ko) 
equation (3.16) is written as 
JO twl2) - gt? dz = FCS), (3.17) 
Following the treatment of Vassell [43], this can be transformed into an Abel integral equation 
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by differentiating with respect to € giving 


HQ as Oe 
de = 8 O=~5 [e(z) — €]/?- 


The independent variable z is eliminated in favor of the dependent variable « by letting 


(3.18) 


f= = [ G(«)dk (3.19) 


where G(x) is an unknown function, and in fact we will not have any reason to determine 
G(«) explicitly. 
Since «(z) is the index of refraction squared, at the reference depth «(0) = 1, and at the 


turning depth «(z;) = &. Equation (3.18) becomes 


n= 5 | (3.20) 


which is an Abel integral equation and can be inverted [28] to give 


: __ 2d f* H(f)a 
Using equation (3.19) and integrating gives the relation 
1 
are | ewe (3.22) 
eel 


which can be used to calculate «(z) from the H(£) obtained by differentiating the function 
F(€) defined by the measured horizontal wavenumbers and associated values of the WKB 
phase integral. The value of the sound speed at one depth is now used to convert the index 


of refraction to a sound speed profile. 


3.2.2 Derivation of inversion method 2 


In the case of an asymmetric sound channel with two turning points, we are able to obtain a 
relation for the distance between the points on the two branches of the sound speed profile 
having the same sound speed. The information contained in the phase integral is insufficient 
to obtain the depths separately. An analagous problem in quantum mechanics is the deter- 


mination of a potential from bound state information and the technique we use is referred 
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to as the JWKB semi- classical approximation in the quantum mechanical context [10],{44]. 


Define the inclusion as 
21 
I(é) =a (x = £)de (3.23) 
z2 
and the excursion as 


AOE = [ 62. — ine (3.24) 


It is the excursion which we will be able to obtain from measurements of the horizontal 
wavenumbers and associated values of the phase integral. For a monotonic profile, the second 
turning point is at the surface (z = 0), and the excursion will give the index of refraction 
versus depth directly, and in asymmetric channels the surface is the upper turning point for 
modes interacting with the surface so that the excursion will be the sound velocity profile for 
a portion of the profile. 


Define the adjoint fractional integral as [39] 
1 : -1 
We" f(a) = ao [ (6-2) FOG; Rev > 0. (3.25) 
(vy) Je 
Applying this to both sides of the WKB phase integral equation with v = —1/2 gives 
Emaz z Emaz 
Pe -9? | PP wae] ae = [OOP Fea. (8.26) 
g Z1 g 
Figure ( 3-1) shows the region over which the square of the index of refraction is being 


integrated. Interchanging the order of integration, it is apparent that the integral over €’ will 


be zero for € greater than x. This leaves 
a We Cis f nue ‘yo . CE) / 
7 dz / (=) ie = | aa int (3.27) 


_ Poe ae 2 [PED ae 
Me) = f° («- az == | eae (3.28) 


Taking the derivative of both sides of this equation with respect to € results in an expression 


ee | PN as | 
x= ff dz = - 5, = | One. (3.29) 


Thus from the values of F(€), we determine the inclusion for the sound speed profile, and 


which reduces to 


for the excursion 


use the inclusion to find the excursion. 


33 


1.00 (| ee as 
inal 


SS i ess 
aS 
Ain hee ee ses 
0.98 o | PS eee) ees 


Kappa 


0.97 


0.95 
| 
0 000 1000 1500 2000 2900 3000 


Depth (m) 


= 


Figure 3-1: The integration over ¢’ is carried out over the cross-hatched region extending 
from € to the curve «(z). 


3.2.3. Examples of analytically invertible profiles 


The n?(z) linear profile and the parabolic profile are two examples of sound speed profiles 
which are analytically invertible with these techniques, and which have well known expressions 
for the eigenvalues. These profiles will serve as test cases for checking the accuracy of the 


intermediate steps in the numerical inversions. 


The n?(z) linear profile 


With the square of the index of refraction given by 
me) = (oh 1 — az. 


the WKB phase integral can be evaluated directly 


F()= f° (az) - gi dz = 20-9). (3.30) 


a 
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Using the first inversion method we take the derivative with respect to € giving 
H() = —A-9?, (3.31) 
The integral expression for depth in terms of the index of refraction is 
2 eal 7 US 
a) eee a d SS V 
z= f (*) (3.32) 


which upon integrating and solving for « returns the original expression for the square of the 


index of refraction. 


The parabolic profile 


For the parabolic profile the square of the index of refraction is given by 


K(z) = n*(z) = 1—- a? 2’, (3.33) 
and the inclusion is 
2ty 4 
1(¢) = f “(= a22?) - gaz = — (1-9). (3.34) 
Zto a 
Taking the derivative with respect to € gives the excursion 
0 2 
> X(€) = —(1- ie (3.35) 


The depth of the turning points is given by the condition 
(1 — a?z”) —- £} =0orz=+(1-£)/?/a 


and the excursion is twice this as expected. 


3.3. Implementation 


In order to determine the sound speed with these methods, we require a set of measurements 
of the horizontal wavenumber in the waveguide and the sound speed for at least one depth ( 
this last permits conversion of the index of refraction profile to a sound speed profile). For 
inversion method 2, we assume that the sound speed at the surface and both the depth an 


value of the minimum sound speed in the profile are known. From the calculations of the WKB 
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phase integral values for the various cases, it is clear that the mode numbers of the measured 
eigenvalues must be correctly identified, and that the boundary interactions experienced by 
each mode must be correctly identified. From the eigenvalues and associated WKB phase 
integral values, the variables £ = (k,,/ko)* and F(€) are computed. One additional point, 
corresponding to horizontally propagating energy, is added at £ = 1 and F( = 0) = 0. 

The first step in either inversion method is to fit the data using either a least squares basis 
spline or a cubic smoothing spline. Either of these two spline methods allows for noise in the 
data. The spline routines essentially provide a piecewise polynomial that provides a best fit 
to the data in a sense which depends on the type of spline routine used. The use of splines is 
convenient in that the integrations and differentiations required for the inversion are easily 
performed once the coefficients for the spline representaion have been computed. All the 
computer routines used for fitting the data with the spline, differentiation, and integration 
are in the IMSL, Inc. MATH/LIBRARY [1] . 

Before discussing the differences in the two spline fitting methods, it is worthwhile to 
define the basis spline which first requires a definition of the divided difference (for a com- 
plete discussion of splines see reference [14]). For a function g which is given at the points 
Zj,---,2j44, the kth divided difference of g at the points z;,...,2;4 1 is the leading coefficient 
of the polynomial of order k +1 which agrees with g at the points 2;,...,2;4,}; it is denoted 
by [z;,...,2i4%]g. The agreement between the two functions referred to here means that if a 
point occurs m times in the sequence z , then the two functions and their derivatives agree 
m-—fold in that the two functions and the derivatives up to the (m—1)th derivative are equal. 


For a sequence (29,21,..-,2N), if we let 
Vice pee — ry) 
and 
kd ne ian (Eeeeia )eos | eS EE Ae ee 7) 


(i.e. the (x — z;) factor is removed before evaluating the product at the point z;), then the 
divided difference is [4] 


N 
(to; meng = Be (3.36) 
j=0 j 
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The basis spline or B-spline is defined as follows [14]: for a nondecreasing sequence t = (t;), 
the ith (normalized) B-spline of order k for the knot sequence ¢ is denoted by B; 4, and is 
defined by the rule 
bet = (tine — te)[ti,-.-, tine]? -— 2g. 

The truncation function (x — t)4 is defined as max{0,z —t} and the -notation is used to 
indicate that z is fixed and (¢ — x), is to be considered as a function of ¢ alone. Note that 
B;x2(2) is 0 for x ¢ [t;, t:4,]. A spline function of order k with knot sequence ¢ is any linear 
combination of B-splines of order k for the knot sequence ¢. 

The least squares spline routine calculates a weighted discrete Lz approximation to the 
given data (2;, f;) for i = 1,2,...,N, (it finds the coefficients a;) to minimize the weighted 
square error between the data and the spline i.e. 

N m 

Ds ie > 45 B;(z:) |? wi. 

i=] j=l 
The number of data points is V, B; is the 7th spline of order k, m is the number of coefficients 
(the number of B-splines making up the polynomial representation), and w; are the data 
weights. The order k is the order of the polynomial pieces ( a polynomial of order k is a 
polynomial of degree & — 1; for a cubic polynomial k = 4). In general, the weights for our 
problem have been selected such that the measured data are equally weighted, and the added 
point at (1,0) is very heavily weighted compared to the measurements. These weights can 
be tailored to fit the confidence in the measurements. 

The smoothing spline is a natural cubic spline approximation with knots at the data 
abscissas where the term ‘natural’ refers to the end condition imposed. In addition to the 
constraints imposed by the required agreement between the function and the spline at the 
knots, some conditions or constraints must be given at the endpoints of the interval on which 
the spline is being calculated. In the case of the natural spline of order k = 2 and N data 


points, the condition specified is 
fO)(a)) = fo Gx) = 0; es ree 
which in the case of a cubic spline (m = 2) becomes 
fC) =4 Gr) =o 
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In our problem, the data error variance is unknown, and the required amount of smoothing 
is determined in the IMSL smoothing spline routine by using the generalized cross validation 
method of Craven and Wahba [13]. Modelling the data as 

y(t) = g(t) + (t);  t € [0,1], 


the problem is to construct g(t) from the discrete data points y;. The function g(t) is assumed 


to be in WS”), where 
wi”) = {g:g” abs. cont., y=0,1,...,m—1, g™eE L,{0, 1}}. 


The estimate of g is gn,,, the solution to the problem: Findf € wi” to minimize 
- > ) -wP +f (Faw) 


The solution gn,, is the a spline of order 2m—1 [38] with the parameter controlling 


the tradeoff between the smoothness of the solution measured by 


l 2 
; / fo) (u)| aa 
0 
and the infidelity to the data measured by 
1 n 
BUGS ae 
j=l 
The generalized cross-validation estimate for the minimization of the average square error is 


defined by [13] : 
V(A) = =I] (FAQ) y IP / [=Trace(I- AQa))) (3.37) 


where y = (y1,.--, Yn)’ and A(A) is the n x n matrix satisfying 


igaexteyn)e cee Soe CEO = A(A)y. 


Once the initial spline fit to the data has been done using either of the two methods, the 
paths of the two inversion methods diverge. The important points related to splines for this 
work are: (1) the spline routines fit piecewise polynomials to the data (the polynomial pieces 
are generally cubic, although the least squares routine allows the order of the polynomial 
pieces to be specified); (2) the spline routines used allow for the presence of noise in the data; 
and (3) the integrations and differentiations are based on the polynomial coefficients from 


the spline fit to the data. 
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3.3.1 Implementation of inversion method 1 


The function H(£) is generated by taking the derivative of the spline representation of F() 
at the data points €; of the input data. A spline is then fitted to H(£). Since noise in the 
original data may be amplified in the processs of computing the derivatives, a smoothing or 
least squares spline is also used at this stage of the computations. The integral that must 
now be evaluated is 


=-5 y es (3.38) 


< €-n)P 


which has a singularity at k. To avoid problems with the singularity during integration, the 


change of variables ({4], [43]) u? = (€ — x) is made and the integral to be evaluated becomes 


(1 x)i/2 
z= -- | H(« + u*) du. (3.39) 
0 


This is evaluated at points over the range of « (which coincides with the range of €) using 
the spline and integration routines. The one required measurement of c(z) is then used to 


convert the points of (z,n(z)) to a sound speed profile c(z). 


3.3.2 Implementation of inversion method 2 


The integral which must be evaluated to obtain the inclusion is 


I(§) = = an aan 7y ae (3.40) 


which is of the same form as equation ( 3.38), and the same change of variables is used to 
eliminate the singularity. The integration is again done using the same spline routines. The 
resulting points of the function J(£) are fitted with a spline and differentiated with respect 


to € to produce the desired result, —X(€). 


3.4 Application of the Inversion Methods to the Analytic 
Profiles 


The two analytically solvable profiles introduced earlier provide excellent examples with which 
to test the inversion methods since the results of the inversion calculations can be compared 


to known values during the intermediate stages of the computations. Expressions for the 
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eigenvalues of these sound speed profiles are well known, and the eigenvalues can be computed 
without the use of a propagation modelling program. In the case of the n?(z) linear profile, we 
take the ocean as a half space with a pressure release surface, and use a = 2.28 x 10°° m7! 
(this gives a sound speed gradient of 0.017sec~! at a depth of 1000 meters corresponding 
to sound speed profile where pressure is the only factor affecting the sound speed i.e. an 
isovelocity, isosaline ocean). For the parabolic profile the ocean is treated as an infinite 


medium, and the coefficient a in the equation for n?(z) is 1.257 x 1074 m7?. 


In each case 
the 50 eigenvalues used have been calculated for a frequency of 220 Hz and a 1500 meter per 
second reference sound speed. Unless otherwise specified, all the calculations have been done 


using the least squares basis splines. 


3.4.1 Method 1 applied to the n*(z) linear profile 


The data points for the function F(£) are created by dividing the input values of the horizontal 
wavenumber by ko = w/co, and squaring the result to give €. The eigenvalues for this example 


were calculated using the expression [8] 
2 = kB — yn (3.41) 


where yn are the Airy function zeros, H = (ak2)~!/3 (equation(2.28)), and a is the coefficient 
in the equation for the index of refraction. The asymptotic expressions [2] were used for the 
Airy function zeros; for the particular parameters used in thes examples this approximation 
did not introduce a significant error. The corresponding value of F'(€) is obtained by dividing 
the assumed value of the WKB phase integral by kg. For this profile, all modes interact with 
the surface and have one turning point, and the WKB phase integral is given by equation 
(3.4) 
[oe - eae = a(n - 1/4) ale), 3 ee 


The point at (1,0) is added to the data set, and the spline coefficients calculated. Figure 
( 3-2) shows the function F() for this particular case. The background curve represents 
F(€) computed using the analytic expression, and the discrete points are the values from 
the inversion calculations showing that there is good agreement between the actual value of 


the phase integral and the assumed value. Using the spline coefficients, the first derivative 
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Figure 3-2: F(£) for the n?(z) linear profile. The discrete points are the input data to the 
inversion program generated from the horizontal wavenumbers and the associated values of 
the WKB phase integral. The curve is F(£) computed using the analytic expression. 


of F(€) with respect to € is computed giving H(€) as shown in figure (3-3). In calculating 
H(€), a noticeable error occurs near £ = 1.0 due to the end conditions in the spline fit for 
F(€) which do not coincide exactly to the behaviour of the function. 

A spline is fitted to the function H(£) represented by the the resulting discrete values of 
H(€), and, using the change of variables u* = (€ — «), the integral 


(1—n)!/2 
z=-= | H(« + u*)du 


w Jo 
is evaluated by varying « through the range for which input data is available. In essence, for 
each chosen value of «k, the depth at which a ray with the local propagation angle defined 
by « is found. As a matter of convenience the range of « or € has been equally divided into 
the number of equal segments corresponding to the desired number of depth values to be 
computed. Figure (3-4) shows the results of the inversion for the n*(z) linear profile using 


50 modes calculated at a frequency of 220 Hz. The convention in presenting sound speed 
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Figure 3-3: H(£) for the n*(z) linear profile. The discrete points are the values of H(£) 
computed by the inversion program. The curve is the function calculated using the analytic 
expression. 


profiles will be to show inversion outputs as discrete points with the true profile plotted as 
a solid curve. The sound speed error as a function of depth for this case is shown in figure 
(3-5). Note that the largest errors occur near the surface where the spacing between the data 
points along the € axis is greatest, and where the extrapolation from the measured data to 
the added point placed at £ = 0 is made. This region also corresponds to the region with the 


greatest error in the calculated values of H(€). 


3.4.2 Method 2 applied to the n?(z) linear profile 


The initial data preparation is the same for both inversion methods since this is a monotonic 
profile. Once the spline has been computed for the function F'(€), the inclusion is computed 


from equation (3.40) using the change of variables u* = ¢’ — € to avoid problems with the 
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Figure 3-4: The sound speed profile for the n*(z) profile recovered using inversion method 1. 
The discrete points are the output values from the inversion. 
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Figure 3-5: Sound speed error as a function of depth for the recovery of the n?(z) linear 
profile using inversion method 1. 


Ad 


singularity at £’ = €. The integral to be evaluated becomes 


K(g) = = = F(u? + du (3.42) 


After computing a spline fit to the points of /(€) generated by the integration, we take minus 
the derivative of /(€) with respect to € to arrive at the excursion as shown in figure ( 3-6). 


The error, plotted as sound speed error as a function of depth, is shown in figure (3-7). 


3.4.3. Method 1 applied to the parabolic profile 


The values of € for the given input data are calculated in the same way as for the previous 
case; however, the values of F(£€) are calculated using equation ( 3.6) 
ko (ey? C27 hel) 2 = od ee 
2to 

since there are two turning points and no boundary interactions for this profile. Because 
the turning depths are equidistant from the channel axis, F(€) is divided by two, and the 
inversion proceeds as before. The resulting sound speed profile is for the region below the 
channel axis, and it is simply a matter of using the profile symmetry to generate the other 
half of the profile. Figure ( 3-8) shows the results of applying inversion method 1 to the 
parabolic profile, and figure ( 3-9) is plot of the inversion error presented as excursion error 
as a function of depth. Excursion error as a function of depth is a clearer way to evaluate 
the inversion results for two branch profiles recovered using method 2, and for consistency 
the errors for inversion method 1 are presented in the same manner. The excursion error in 
the case of a symmetric profile recovered using inversion method 1 is twice the error of the 


output for the single branch of the profile. 


3.4.4 Method 2 applied to the parabolic profile 


The parabolic profile involves modal interactions with two turning points, and inversion 
method 2 gives only excursion as a function of sound speed. Given prior knowledge that 
the profile is symmetric this may, if desired, be converted to a sound speed profile using 
symmetry as in inversion method 1. Figure (3-10) is the excursion for the parabolic profile 
computed using inversion method 2, and figure ( 3-11) plots excursion error as a function of 


sound speed for this inversion. 
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Figure 3-6: Excursion (depth) versus sound speed for the n*(z) linear profile calculated using 
inversion method 2. 
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Figure 3-7: Sound speed error versus depth for inversion method 2 applied to the n?(z) linear 
profile. 
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Figure 3-8: Depth versus sound speed for the parabolic profile recovered using inversion 


method 1. The points below the channel axis are the actual inversion results, and the profile 
symmetry was used to generate the profile above the axis. 
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Figure 3-9: Excursion error as a function of sound speed for inversion method 1 applied to 
the parabolic profile. Excursion error is twice the depth error. 


49 


200 


400 


600 


800 


m) 


1000 


Excursion ( 


1200 


1400 


1600 


1800 


1498 1500 1502 1504 1506 1508 1510 


Sound speed (m/s) 


Figure 3-10: The excursion for the parabolic profile computed using inversion method 2. The 
excursion is the distance between points of equal sound speed on the two branches of the 


sound speed profile. 
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Figure 3-11: Excursion error versus sound speed for inversion method 2 applied to the 
parabolic profile. 
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3.5 Application of the inversion methods to a bilinear profile 


The sound speed profile in the ocean is determined by the temperature, pressure, and salinity 
which in general leads to sound speed profiles more complicated than the ideal profiles we have 
considered. A nominal deep ocean sound speed profile consists of a region near the surface, 
the thermocline, where the temperature effects dominate and the sound speed decreases from 
its surface value to a minimum at the axis of the deep sound channel. Below the depth of 
minimum sound speed, the increase of pressure with depth is the dominant influence, and 
the sound speed increases with depth. Although the sound speed structure can be quite 
complicated, we will use a simple bilinear model with a channel axis at 686 meters (figure 
(3-12) to further illustrate the inversion methods. 

The modal eigenvalues for this model were computed for a frequency of 220 Hz using the 
SACLANTCEN normal-mode acoustic propagation model program (SNAP)(27]. Initially 
the data set was truncated before applying inversion method 2 in order to use only modes 
whose upper turning point is deeper than 100 meters (to avoid effects of interactions with the 
surface). The resulting excursion as a function of sound speed is shown in figure ( 3-13), and 
the excursion error is shown in figure (3-14). The greatest error occurs near the channel axis 
where the profile smoothness assumptions are not met, and where the extrapolation from the 
measured data to the assumed point at € = 1.0 is made. 

Next we use a data set containing the eigenvalues for normal modes with a lower turning 
depth shallower than 3200 meters. The initial portion of this data set consists of eigenvalues 
for modes with two turning points while the eigenvalues at the end of the data set are from 
modes that interact with the surface and have a lower turning point in the water column. 
We assume that the value of the sound speed at the surface is available which allows us 
to determine the value of the horizontal wavenumber separating the two types of modes by 


setting the WKB phase integral integrand equal to zero so that 


= 


ae (it) 
For the modes with two turning points F'(€) is calculated using equation (3.6) 


F(€) = 1(n — 1/2)/ko, 
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Figure 3-12: The bilinear sound speed profile model. 
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Figure 3-13: The excursion for the bilinear profile using eigenvalues for modes with an upper 
turning point deeper than 100 meters. 
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Figure 3-14: Excursion error for the bilinear profile in the region with two turning points. 


and for the modes with only one turning point equation (3.4) gives 
F(§) = x(n — 1/4)/ko. 


Applying inversion method 2 gives the excursion versus sound speed which in the portion 
of the profile having two turning points is the difference in depth between points of equal 
sound speed, and in the portion of the profile with a single turning point is the depth at 
which the value of sound speed occurs. Figure (3-15) shows the excursion versus sound speed 
obtained. As before, large errors occur near the channel axis, plus there are large errors 
at the depth where the transition from two to one turning point is made. This latter error 
appears to result from errors in the spline fit at the point where the form of F'(€) changes. If 
we include a sediment layer with a 1 m/s sound speed gradient starting at a depth of 4000 
meters, then the inversion using method 2 will give the distance between the two turning 
points of equal sound speed for sound speeds in the range from 1478 m/s to 1507 m/s and 
the actual depth at which the sound speed occurs for sound speeds greater than 1507 m/s 


(figure ( 3-16)). Consequently we obtain the sediment sound speed profile directly and over 


00 


a portion of the water column we also obtain the sound speed profile directly. Figure (3-17) 
shows the excursion error over the entire profile. The regions of largest error occur where the 
form of either the sound speed profile or the function F(é) changes. The outputs presented 
here assume that in addition to the modal eigenvalues, the sound speed at the surface and 
the minimum sound speed in the profile are available. If the last item is not known, it may be 
possible to use an iterative procedure to find an approximation to the excursion as a function 
of sound speed (i.e. find the minimum sound speed) or to separate the portions of the sound 
speed profiles above and below the channel axis; however, we have not pursued this. The 
final approach for this profile is to treat it as a symmetric profile by defining an ‘equivalent 


symmetric profile’ [35] such that 

resp = 5 (/24 [+1 201) (3.44) 
where z, and z_ are the distances above and below the channel axis to points having sound 
speed c. Inversion method 1 is applied to the data for the region with two turning points 
(F(€) must be divided by 2 as in the case of the parabolic profile) to give the portion of the 
equivalent symmetric profile below the channel axis, and symmetry is used to generate the 
portion above the axis. Figure ( 3-18) illustrates the results of using this method which, in 
a sense, is averaging the vertical phase contributions of the two portions of the sound speed 


profile. 
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Figure 3-15: Excursion versus depth for the bilinear profile using eigenvalues for modes 
turning above 3200 meters. Below 2813 meters this curve is the sound speed profile. 
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profile. 
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3.6 Performance of the inversion methods 

For any inverse method, questions of practical interest concerning the solution are: 
e Is the solution unique ? 
e What is the resolution of the inversion? 
e How well does the inversion method perform with imperfect or noisy data? 


3.6.1 Uniqueness 


The question of uniqueness for Sturm-Liouville problems has been extensively studied ( see for 
example [29], [32]), and the requirements for unique recovery of the coefficients are well known. 
If we model the problem as a horizontally stratified medium consisting of the water column, 
fluid sediment layers, and a semi-infinite halfspace which can be considered as providing an 


impedance boundary condition, then the problem can be written on a finite interval [0,1] as 


d?y 

as? + (A - q(s))v =0 (3.40) 
sin av(0) + cosav’(0) = 0 (3.46) 
sin Bv(1) + cos Bu'(1) = 0 (3.47) 


The general problem is to find a,Z, and q from spectral data. The basic data set is 
the infinite set of eigenvalues A; < A2 < ... for the given Sturm-Liouville problem. If the 
potential is symmetric about s = 1/2 and a = 7 — @ this is sufficient to provide for unique 
recovery of a and q; however, for a general profile, an additional infinite sequence of data will 


be required. Three different possibilities for this second set of data include [32] 


1. a second set of eigenvalues M < Ao <.... where a different boundary condition is 


applied in place of equation (3.47), i.e. 8 would be replaced by B: 


2. the set of normalization constants 6, =|| vp ||? /[vn(0)]?,n =1,2,..., whenO0<a<zr 


or by =|| vn ||* /[vi,(0)]? when a = 0. vp is the eigenfunction for the eigenvalue A,,, and 


| on II? = Jo [on(s)]?ds; 
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3. the set of ratios vz(1)/v,(0) when 0 < a,f < = or v),(1)/v,(0) when a = 2 = 0, or 


similar ratios fora =0,0< B<xzorB=0,0<a<r. 


In addition, Hochstadt and Lieberman [26] show that uniqueness is possible given the set of 
eigenvalues \j < A2 < ... and knowledge of q(s) on the interval [1/2, 1]. 

For the problem considered here, we have one partial set of eigenvalues and knowledge of 
one boundary condition (the pressure release surface). Conseqently, we cannot expect that 
the answer provided by the inverse methods will be unique. From a practical standpoint this 
means that there may be additional information on a finer verticle scale than the frequency of 
our experiment could measure, but which could be obtained at a different frequency; or there 
may be a variety of models possible for the region below the depth of deepest penetration 
of the rays equivalent to the modes used for the inversion which would give the same finite 
set of eigenvalues as we have measured. It should be noted that, if horizontal wavenumbers 
for two (or more) frequencies are used, the transformation to the variable € will map all the 
eigenvalue data to the inerval [0,1], and all the data will fit on the same curve F(£) since 


F(€) is normalized. This allows us to combine data obtained using two frequencies. 


3.6.2 Resolution 


One of the advantages of perturbative inversion methods is that the linear inverse theory 
used allows a quantitative estimate of the resolution provided by the data (for a discussion of 
linear inverse theory see [33]). With these nonlinear inversion methods we can only provide 
a qualitative description of the depth resolution. Inversion method 1 uses the function H(€) 
obtained by interpolating and differentiating the input data to find the depth at which a 
particular value of the index of refraction will cause a ray with angle sin~!./€ to turn. We 
can expect that the resolution will be related to the depth difference between turning depths 
for successive modes. Better resolution should be obtained if the sound speed profile and 
the frequency are such that the turning depths of adjacent modes are closely spaced. The 
turning depths are found using the condition 


WwW 


and high frequencies or large sound speed gradients give more closely spaced turning depths. 


In integrating over « to obtain the depth for a given index of refraction, we divided the range 
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of « uniformly, and the output values from the inversion will be closer together in portions 
of the sound speed profile with the highest gradients. This can be seen in figure (3-16). 
Inversion method 2 approaches the problem in a similar manner finding the difference 
in depth between points on the profile having the same sound sae and similar arguments 
concerning the spacing of the mode turning points pertain to the depth resolution of this 


method. 


3.6.3 Performance of the inversion methods in the presence of noise 


Because of the nonlinearity of the problem, it is difficult to obtain an analytic expression 
giving an indication of how the inversions will perform when the input data is not exact. For 
inversion method 1, we obtained (equation (3.18)) the Abel integral equation 


GAS eo ee ey dz 
eM) 5h [x(z) — §]!/? 


and a similar Abel equation was obtained in method 2 where, rather than inverting the 
equation, we differentiated with respect to € to obtain the desired result, equation ( 3.29) 
(the excursion). We have considered that the values of F'(€) are exact, and, while some noise 
will be introduced in differentiation due to the inexactness of the values € associated with 
F(€), the noise enters primarily in the kernel of the integral equation. Thus we will approach 
the question of noise effects by numerically testing the inversions using the n?(z) linear profile 
as a test case. 

In carrying out an experiment to measure the eigenvalues in the ocean waveguide, the 
environment and the instruments will introduce noise into the data, the process of computing 
the eigenvalues will have some inherent resolution limit imposed by array size and spatial 
sampling factors, and the computation process will introduce numeric inaccuracies such as 
round-off error. Rather than attempt to quantify all these factors which will vary from 
experiment to experiment, we will assume that we have data with random errors of known 
maximum value. This noise will be modelled as an additive roundoff error with a uniform 
probability distribution from —5.0 x 107” to 5.0 x 10~” where n is the largest decimal position 
in the eigenvalue affected by the noise. The two inversion methods were run with varying 
error magnitudes using both the least squares basis spline routines and the cubic smoothing 


spline routines relying on cross validation to choose the smoothing parameter. Figures (3-19) 
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Figure 3-19: RMS error for inversion method 1 using least squares basis splines with noisy 
data. 


through figure (3-26) plot the root mean square error and the absolute value of the maximum 
inversion errror as functions of the logarithm of maximum introduced error. The results show 
that the smoothing spline can handle a greater magnitude of noise, and that both inversion 
methods can handle reasonable amounts of noise in the data. If a mode is missed, 
i.e. the mode spacing and noise are such that one mode is missed and all subsequent modes 
are incorrectly numbered by one, then the errors in F(€) will produce errors in the output. 
Figure ( 3-27) illustrates the error produced when the eigenvalue for mode 25 of the n*(z) 
linear profile is deleted, and the remaining 24 eigenvalues misnumbered. The sudden jump 
in F'(€) leads to a sharp jump in the error of the output at the depth corresponding to the 


missing mode, and then the error decreases tending to a constant value. 
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Figure 3-20: Absolute value of the maximum error for inversion method 1 using least squares 
basis splines with noisy data. 
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Figure 3-21: RMS error for inversion method 2 using least squares basis splines with noisy 
data. 
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Figure 3-22: Absolute value of the maximum error for inversion method 2 using least squares 
basis splines with noisy data. 
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Figure 3-23: RMS error for inversion method 1 using smoothing splines with noisy data. 
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Figure 3-24: Absolute value of the maximum error for inversion method 1 using smoothing 
splines with noisy data. 
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Figure 3-25: RMS error for inversion method 2 using smoothing splines with noisy data. 
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Figure 3-26: Absolute value of the maximum error for inversion method 2 using smoothing 
splines with noisy data. 
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Figure 3-27: Sound speed error for the n?(z) linear profile with mode 25 eigenvalue removed, 
and higher modes incorrectly identified. 
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Chapter 4 


Application to a Shallow Water 
Waveguide 


In order to illustrate the application of this inversion technique to a specific problem in 
underwater acoustics, we will consider the problem of determining the acoustic properties 
of the sediments in a shallow water waveguide. The objective is to provide a portion of the 
geoacoustic model for the waveguide. A complete geoacoustic model [25] includes information 
on the geology and topography of the bottom, density profiles, shear wave properties, and 
compressional wave properties. This inverse method is capable of providing estimates for 
compressional wave speed and density profiles, although we can only illustrate application to 
finding compressional wave speed profiles. | 

The data for the inversion is obtained from an experiment [20], [30] in which a continuous 
wave source at a fixed depth is towed away from a pair of moored receivers which record 
sound field pressure (amplitude and phase) as a function of range (see figure 4.1). In actually 
carrying out the experiment, one source operating at two different frequencies is used which, 
in principle, allows determination of the density as well as sound speed. The pressure field is 
sampled at least every half wavelength in range over an aperture of several kilometers, and 
estimates of the horizontal wavenumbers for the propagating sound field are obtained from 
the pressure field measurements either by using Prony’s method [15] which models the field as 
a sum of complex exponentials, or by numerically Hankel transforming the measurements to 
obtain the depth dependent Green’s function. In previous work, the bottom properties have 
been determined using the data from this type of experiment [21] by assuming a geoacoustic 


model, computing the theoretical Green’s function for the model, and calculating the mean 
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Figure 4-1: Experimental configuration for measuring the eigenvalues in a shallow water 
waveguide. The parameters shown are for an experiment in Nantucket Sound, Massachusetts 
[21] 


square difference between the two Green’s functions. The model is varied and the process 
repeated until the total mean square difference over the four frequency /receiver depth combi- 
nations is a minimum. Rajan et al. [37] apply perturbation techniques to the problem using 
the eigenvalues obtained via the Hankel transform of the pressure field as the input. 

In constructing the geoacoustic model with data from this experiment, information on 
water column properties will be available from separate measurements to provide the sound 


speed value needed to transform the index of refraction output into a sound speed profile. 


4.1 Bottom Models and Inversions with Synthetic Data 


We will use two model waveguides each having a one hundred meter deep isovelocity water 
column, a one hundred meter thick sediment layer, and a basalt basement (half-space) having 


density, compressional wave speed, and attenuation equal to nominal values for basalt (shear 
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properties are not included in the models). Two models will be used for the sediment layer: 
a model with a terrigenous sediment layer made up of silty sand, and model consisting of fine 
sand sediment. The data on which both models are based is contained in Hamilton’s 1980 
Teview article on geoacoustic models [25]. For both bottom types the sediment density has 
been taken as constant due to constraints in SNAP, the acoustic modelling program which 


was used to calculate the horizontal wavenumbers for the sound field in the waveguides [27]. 


4.1.1 The Terrigenous Bottom Model 


The constituents of terrigenous sediments are derived from weathering and erosion of rocks 
found mostly on land [40]. Hamilton obtained the following regression equation for velocity 
as a function of the depth in the sediments using data assembled for twenty areas in upper, 


unlithified layers of mostly turbidites in various oceans [25] 
V(z) = 1.502-+41.301z — 0.7412? +.0.2572°. (4.1) 


The depth z is in kilometers, and the velocity V is in kilometers per second. The density in 
the sediment layers has been taken as p = 1.0g/cm? so that the eigenvalues will be governed 
solely by the sound speed profile and the boundary conditions at the ocean surface and the 
sediment/basement interface. Water column sound speed has been set equal to the sediment 
sound speed at the water-sediment interface for the same reason. Figure ( 4-2) is the sound 
speed profile for the waveguide with a terrigenous bottom model. 

SNAP produced eigenvalues for 53 propagating modes in this waveguide; however, only 
the first 15 of these eigenvalues were used in the inversion in order to avoid effects arising 
from interactions of the modes with the basement half-space since we are interested in testing 
our ability to recover the sediment sound speed profile. The highest mode used turned at a 
depth of 176 meters. The modes were treated as one turning point modes, and both inversion 
methods were used. Figures ( 4-3) through ( 4-6) show the inversion results with the errors 
plotted as depth error versus sound speed. As expected, the largest errors occur at the 
wWater-sediment interface where the extrapolation from the data point corresponding to the 
first mode to the point at € = 1, F(€) = 0 is made and where the profile does not meet the 
smoothness assumptions. We have also ignored the possibility of reflections at the interface 


which may introduce errors in F(£) for the low order modes (those with large equivalent angles 


v1 


of incidence). The apparent increase in error at the lower depths reached by the inversion 
is believed to result from inaccuracies in the spline fit near the end of the data. Note that 
the inversion does not produce any points on the isovelocity portion of the profile since there 
are no turning points in an isovelocity segment (the one inversion point shown at the water 
surface is the reference sound speed i. e. the assumed minimum in the profile). Terrigenous 
sediments are predicted to have a sound speed at the water/sediment interface lower than 
the water sound speed [25]. The inversion methods cannot distinguish such a low velocity 
zone because they are finding turning depths which result from increasing sound speeds. In 
figure (4-7) are the results for model 1 with the water column sound speed increased to 
1545 meters per second to give a low velocity zone. The reference sound speed used in the 
inversion was taken as 1545 meters per second which, because of the normalization of the 
eigenvalues, produced a value of € > 1 for the first eigenvalue and, to produce the results of 
figure (4-7), the first eigenvalue was discarded. Thus the data give a clear indication that a 
grossly inaccurate reference sound speed has been used. If the first eigenvalue is retained and 
the correct reference sound speed of 1511 meters per second used, then the result obtained is 
shown in figure ( 4-8). These results indicate that the inverse methods will at least indicate 
the presence of a low velocity zone at the water/sediment interface, and it may be possible 
to gain some information about the value of the minimum sound speed through iterative 


adjustments to the reference sound speed. 
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Figure 4-2: The shallow water waveguide model with a terrigenous sediment bottom. 
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Figure 4-3: Inversion results for waveguide model 1 using inversion method 1. The back- 
ground curve is the actual sound speed profile in the waveguide. 
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Figure 4-4: Depth error for inversion method 1 applied to waveguide model 1. 
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Figure 4-5: Inversion results for waveguide model 1 using inversion method 2. The back- 
ground curve is the actual sound speed profile in the waveguide. 
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Figure 4-6: Depth error for inversion method 2 applied to waveguide model 1. 
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Figure 4-7: Inversion results for a terrigenous waveguide with a low velocity zone at the 
water /sediment interface. The eigenvalue of the lowest mode was discarded since it produced 
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Figure 4-8: Inversion results for a terrigenous waveguide with a low velocity zone at the 
water/sediment interface obtained using the correct reference sound speed of 1511 meters 


per second. 
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4.1.2 The Fine Sand Bottom 


The fine sand bottom has a large discontinuity in the compressional wave speed at the 
water/sediment interface and large sound speed gradients in the top portion of the sediment 
layer. A regression equation given by Hamilton for the sound speed in the first twenty meters 


of a fine sand bottom is 


V(z) = 1806.02°°°. (4.2) 


The depth z is in meters and the resulting sound speed V is in meters per second. For 
depths greater than twenty meters, the profile has been extended as a straight line with slope 
1.4168 s—' which is the sound speed gradient of equation (4.2) at a depth of twenty meters. 
Figure ( 4-9) is the sound speed profile for the fine sand bottom model. 

For this model SNAP generated eigenvalues for 47 propagating modes of which 23 were 
used for the inversions (the turning depth of the highest mode used was 174 meters). Both 
inversion methods were again used for the inversions, and, in addition, both the least squares 
and the smoothing spline routines were used. The results are shown in figures (4-10) through 
(4-17). An attempt was made to treat the modes with equivalent incidence angles greater 
than the critical angle as interacting with a critically reflecting bottom which was modeled 
as a 1727 meter per second half-space (i.e. they were treated as case 3 modes); however, 
this made the agreement between the inversion results and the actual profile worse than for 
previous inversions which treated all the modes as simple one turning point modes (figure 
( 4-18)). 

Finally, we used all the eigenvalues in the inversion process for this model (figure ( 4-19)). 
The inversion output sound speed approaches the basement sound speed value while the 


depth value does not depart significantly from the sediment /basement interface depth. 
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Figure 4-9: The shallow water waveguide model with a fine sand bottom. 
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Figure 4-10: Inversion results for waveguide model 2 using inversion method 1 (least squares 
splines). 


81 


4.00 


2.00 
el 0.00 
= 
5 -2.00 
i= 
eb) 
4.00 
Sh 
eb) 
a 

—6.00 

—8.00 

1500 1600 1700 1800 1900 2000 


Sound speed (m/s) 


Figure 4-11: Depth error for inversion method 1 (least squares splines) applied to waveguide 
model 2. 


82 


80 


~ * SOE ieee 
etek x fie Be 


100 


28 


Depth (m) 


140 


160 


180 





1500 1600 1700 1800 1900 2000 


Sound speed (m/s) 


Figure 4-12: Inversion results for waveguide model 2 using inversion method 1 (smoothing 
splines). 
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Figure 4-13: Depth error for inversion method 1 (smoothing splines) applied to waveguide 
model 2. 
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Figure 4-14: Inversion results for waveguide model 2 using inversion method 2 (least squares 
splines). 
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Figure 4-15: Depth error for inversion method 2 (least squares splines) applied to waveguide 
model 2. 
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Figure 4-16: Inversion results for waveguide model 2 using inversion method 2 (smoothing 
splines). 
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Figure 4-17: Depth error for inversion method 2 (smoothing splines) applied to waveguide 
model 2. 
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Figure 4-18: Inversion results for waveguide model 2 using inversion method 1 (least squares 


splines) treating the low order modes as reflecting from a half-space with a 1727 meter per 
second sound speed. 
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Figure 4-19: Inversion results for waveguide model 2 using inversion method 1 (least squares 
splines) using the eigenvalues of all propagating modes. 
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4.2 Experimental Results 


An experiment was performed in the Gulf of Mexico near Corpus Christi, Texas during 
September, 1985 using a configuration similar to that of Figure (4.1). For the Corpus Christi 
site the water depth was 30.3 meters, source depth was 22.86 meters, and the two receivers 
were 1.5 meters and 15.0 meters above the bottom (G. V. Frisk pers. com. 1988) Geologic 
data {31] for the region predicts a 17.5 meter thick layer of silty clay with an estimated sound 
speed at the water/sediment interface of 1517.6 meters per second increasing to 1531.5 meters 
per second at the layer bottom. Since the sound speed in the water column was 1545 meters 
per second, this gives a low velocity zone. Beneath the silty clay layer, a layer of very fine 
sand was predicted to extend to 45 meters with a sound speed of 1737.1 meters per second at 
the interface between the layers and increasing to 1763.0 meters per second at the bottom of 
the layer. In the absence of better information, parameters for this segment were extended to 
a depth of 100 meters. It should be emphasized that these values are based on a study of the 
geology of the region, and the sound speed predictions are based on the type of sediments and 
the estimates of Hamilton’s work for such sediment types [25]. Figure ( 4-20) is the predicted 
sound speed profile in the sediments. The eigenvalues were obtained by S. D. Rajan (pers. 
com. 1986) from the depth dependent Green’s function of the measured sound field via the 


numerical Hankel transform, and are presented in Table (4.1). Water column sound speed at 


_Pieyueney (la) eds wales Mayeoveiie (ie 
90.0 1 0.1972886 
wz 0.182872 


0.5678598 
0.5572134 
0.5440619 
0.5252742 
0.5102440 


140.05 


ao fk GNM 


Table 4.1: Experimentally measured eigenvalues for the Corpus Christi site. 


the experiment site was 1545 meters per second. Since a total of only seven eigenvalues at the 


two frequencies is available, we combined the data to form a single data set which is possible 
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Figure 4-20: Sound speed profile for the Corpus Christi area based on geologic information. 
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Figure 4-21: Corpus Christi area sediment sound speed profile obtained using all seven eigen- 
values at 50 Hz and 140.05 Hz. 


since F(£€) as we have defined it is frequency independent. The results of the inversion using 
method 1 with an inversion program using the least squares spline routines of the PORT 
Library [18] are shown in figure (4-21). (The PORT Library program was used because it 
seemed to handle the small data set better than the program using the IMSL routines.) The 
large oscillation at the bottom of the inversion results suggets a data error at the end of the 
data set (smallest £). After € is calculated, the second 50 Hz eigenvalue is the last data point. 
and if this eigenvalue is removed, the oscillation is also removed as shown in figure ( 4-22)) 
which also shows the results of a perturbative inversion done by S$. D. Rajan (pers. com. 


1986), and the predicted sound speed profile. 
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Figure 4-22: Corpus Christi sediment sound speed profile obtained with inversion method 


1 after removing the second 50 Hz eigenvalue, the profile obtained perturbatively by Rajan 
(solid curve), and predicted profile (dotted curve). 
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Chapter 5 


Conclusions 


5.1 Summary 


In this thesis we have shown that using a WKB inversion technique based on the WKB 
phase integral equations and using modal eigenvalues as input data provides a method for 
obtaining information about index of refraction profiles in the ocean waveguide.Information 
can be obtained for the water column and for the sediment layers. Two approaches for 
transforming the WKB phase integral to a relation between depth and the index of refraction 
were developed. Both methods assume that the value of the WKB phase integral can be 
accurately estimated based on the number of turning points for a mode and/or the types of 
boundary interactions the mode experiences, and both methods rely on knowledge of at least 
the speed of sound at the surface. 

Inversion method 1 differentiates the WKB phase integral equation to obtain an Abel 
integral equation which is then inverted to provide an equation for depth [43] and is only 
usable for monotonic or symmetric sound speed profiles. In the case of a symmetric profile 
the knowledge that the phase contributions above and below the axis are equal allow the 
original WKB phase integral to be divided into two equal portions which essentially puts the 
problem in monotonic form. With a value of sound speed at one depth, the monotonicity 
or symmetry allows conversion of the depth/index of refraction relation to a sound speed 
profile. For a general asymmetric profile, it is possible to to treat the data as though it was 
obtained from a symmetric profile and to obtain an ‘equivalent symmetric profile’ [35]. The 
utility of this approach remains to be seen. 


Inversion method 2 uses fractional calculus methods [39] to find the inclusion (equation 
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( 3.28)) from the WKB phase integral. The inclusion is differentiated to obtain the depth 
difference between two points with the same index of refraction (the excursion). For a mono- 
tonic profile this gives the depth directly and for a symmetric profile the phase integral can 
be split into two equal parts as it was for inversion method 1 allowing determination of the 
full profile. For an asymmetric profile this method gives different answers in the portions 
of the profile with one and two turning points. In the portion with two turning points, the 
inversion result is the distance between two points of equal index of refraction, and in the 
single turning point case the result is the index of refraction. This assumes that we have at 
our disposal the sound speed at the surface (to allow the modes to be separated into the one 
and two turning point groups) and the minimum sound speed in the profile. The physical 
interpretation of the variable € as sin*@ which constrains it to the range (0, 1], suggests that 
it may also be possible to approximate the minimum sound speed value based on the range 
of the normalized eigenvalues. Since a known portion of the inversion result for a profile 
with both one and two turning point modes is the index of refraction as a function of depth, 
it may be possible to use an iterative procedure to arrive at a reasonable excursion versus 
sound speed curve and, perhaps, to separate the two branches of the sound velocity profile 
in the two turning point section, that is to find a good approximation to the entire sound 
speed profile. 

Both methods have been implemented using splines for differentiation and integration. 
The techniques were applied to three idealized profiles to demonstrate the ability of the 
inversion methods to recover the profiles. A significant disadvantage of these methods is the 
present lack of quantitative statements concerning depth resolution and performance with 
inaccurate data. Performance of the methods with inaccurate data was demonstrated with 
numerical experiments using the n*(z) profile after adding random noise to the eigenvalues. 
The inversions are capable of giving results even with some noise in the data, and the use of 
inversion method 1 on a set of experimental data was demonstrated in Chapter 4. The results 
of the inversion of the experimental data show that the inversion results may give a fairly 
clear indication ( large depth oscillations in the output) when an eigenvalue is incorrect. The 
depth resolution is related to the depth separation between succesive turning points; however, 


we do not yet have a quantitative expression for the resolution. 
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In Chapter 4, both methods were applied to two models for a shallow water waveguide to 
illustrate the capability to provide estimates of sediment compressional wave speed profiles in 
shallow water. Finally, inversion method 1 was applied to a small data set from an experiment 
performed in the Gulf of Mexico. Only seven eigenvalues at two frequencies were available, 
and it appears that one value is not accurate based depth oscillations at the end of the sound 
speed range. Although our result is much different from the profile predicted on the basis of 
the geology in the area, it is in general agreement with a result obtained perturbatively by 


S. D. Rajan (pers. com. 1986). 


5.2 Future work 


A number of issues related to this work need to be addressed more fully in the future. We have 
largely ignored the discontinuity in sound speed at the water/sediment interface. In Chapter 
4 we tried treating the low order (large angle of incidence) modes as being critically reflected 
at the interface, but this made the agreement between the inversion and the true profile worse 
than ignoring the presence of the discontinuity. The manner in which discontinuities in the 
index of refraction profile should be treated in the context of WKB theory, and the effect of 
the discontinuities on the cumulative phase needs further investigation. 

In Chapter 4 we alluded to use of the inversions to estimate the density profile as well as 
the compressional wave speed profile. If the effects of density (assumed to be a function of 
depth only) are included, the homogeneous Helmholtz equation becomes [8] 

V? p(n, 2) — Volz) V(r 2) + &(2)a(r 2) = 0, (5.1) 


and separation of variables gives for the depth equation 


Substituting v(z) = P(z)//p(z) gives the Schrodinger type equation [36] 





Fv(2) + (k?(z) + p(z) - a v(z) = 0 (5.3) 
where / 
Bn sil (ead 
ACNE p = Foorrae. (5.4) 


Ht 


The WKB phase integral is now 
ko [ (n2(2) ah aes ¢) dz. 
0 Ko 

An open question is whether we can obtain a depth profile for a new variable y(z) = n*(z) + 
u(z)/ko at two frequencies and then eliminate n?(z) from the pair of equations since the index 
of refraction is frequency independent. This would leave an equation for (z) which through 
equation (5.4) provides a differential equation for p(z). The numerical modelling program 
used in this work required constant densities in the sediments and the basement and did not 
provide a means to pursue the density question. 

The inverse method was applied to one data set which had only six good eigenvalues. 
Identification of existing experimental data from other experiments obtained with arrays 
capable of resolving a larger number of modes (at least 10 to 15) would provide an opportunity 
to more fully test the inversions with experimental data. If datais available at two frequencies, 
use of Prony’s method [15] to find the eigenvalues will provide an estimate of the modal 
attenution in adition to the eigenvalues. Using these inversion methods in conjunction with 
Prony’s method will allow a geoacoustic model incorporating sediment sound speed, density, 


and attenuation to be developed. 
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